{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 20\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dropping pennies\n",
    "\n",
    "I'll start by getting the units we need from Pint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "second"
      ],
      "text/latex": [
       "$\\mathrm{second}$"
      ],
      "text/plain": [
       "<Unit('second')>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = UNITS.meter\n",
    "s = UNITS.second"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And defining the initial state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>y</th>\n",
       "      <td>381 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v</th>\n",
       "      <td>0.0 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "y             381 meter\n",
       "v    0.0 meter / second\n",
       "dtype: object"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "init = State(y=381 * m, \n",
    "             v=0 * m/s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Acceleration due to gravity is about 9.8 m / s$^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "9.8 meter/second<sup>2</sup>"
      ],
      "text/latex": [
       "$9.8\\ \\frac{\\mathrm{meter}}{\\mathrm{second}^{2}}$"
      ],
      "text/plain": [
       "9.8 <Unit('meter / second ** 2')>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = 9.8 * m/s**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'll start with a duration of 10 seconds and step size 0.1 second."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "10 second"
      ],
      "text/latex": [
       "$10\\ \\mathrm{second}$"
      ],
      "text/plain": [
       "10 <Unit('second')>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_end = 10 * s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.1 second"
      ],
      "text/latex": [
       "$0.1\\ \\mathrm{second}$"
      ],
      "text/plain": [
       "0.1 <Unit('second')>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dt = 0.1 * s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we make a `System` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>y             381 meter\n",
       "v    0.0 meter / secon...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>g</th>\n",
       "      <td>9.8 meter / second ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>10 second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>dt</th>\n",
       "      <td>0.1 second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "init     y             381 meter\n",
       "v    0.0 meter / secon...\n",
       "g                                  9.8 meter / second ** 2\n",
       "t_end                                            10 second\n",
       "dt                                              0.1 second\n",
       "dtype: object"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = System(init=init, g=g, t_end=t_end, dt=dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And define the slope function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func(state, t, system):\n",
    "    \"\"\"Compute derivatives of the state.\n",
    "    \n",
    "    state: position, velocity\n",
    "    t: time\n",
    "    system: System object containing `g`\n",
    "    \n",
    "    returns: derivatives of y and v\n",
    "    \"\"\"\n",
    "    y, v = state\n",
    "    g = system.g    \n",
    "\n",
    "    dydt = v\n",
    "    dvdt = -g\n",
    "    \n",
    "    return dydt, dvdt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's always a good idea to test the slope function with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 meter / second\n",
      "-9.8 meter / second ** 2\n"
     ]
    }
   ],
   "source": [
    "dydt, dvdt = slope_func(system.init, 0, system)\n",
    "print(dydt)\n",
    "print(dvdt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're ready to call `run_ode_solver`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>The solver successfully reached the end of the...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                                                 True\n",
       "message    The solver successfully reached the end of the...\n",
       "dtype: object"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results, details = run_ode_solver(system, slope_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>v</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.0</th>\n",
       "      <td>381 meter</td>\n",
       "      <td>0.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.1</th>\n",
       "      <td>380.951 meter</td>\n",
       "      <td>-0.9800000000000001 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.2</th>\n",
       "      <td>380.80400000000003 meter</td>\n",
       "      <td>-1.9600000000000002 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.3</th>\n",
       "      <td>380.559 meter</td>\n",
       "      <td>-2.9400000000000004 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.4</th>\n",
       "      <td>380.216 meter</td>\n",
       "      <td>-3.9200000000000004 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                            y                                   v\n",
       "0.0                 381 meter                  0.0 meter / second\n",
       "0.1             380.951 meter  -0.9800000000000001 meter / second\n",
       "0.2  380.80400000000003 meter  -1.9600000000000002 meter / second\n",
       "0.3             380.559 meter  -2.9400000000000004 meter / second\n",
       "0.4             380.216 meter  -3.9200000000000004 meter / second"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>v</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>9.6</th>\n",
       "      <td>-70.58399999999979 meter</td>\n",
       "      <td>-94.08000000000003 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9.7</th>\n",
       "      <td>-80.0409999999998 meter</td>\n",
       "      <td>-95.06000000000003 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9.8</th>\n",
       "      <td>-89.5959999999998 meter</td>\n",
       "      <td>-96.04000000000003 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9.9</th>\n",
       "      <td>-99.24899999999981 meter</td>\n",
       "      <td>-97.02000000000004 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10.0</th>\n",
       "      <td>-108.99999999999982 meter</td>\n",
       "      <td>-98.00000000000004 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                              y                                  v\n",
       "9.6    -70.58399999999979 meter  -94.08000000000003 meter / second\n",
       "9.7     -80.0409999999998 meter  -95.06000000000003 meter / second\n",
       "9.8     -89.5959999999998 meter  -96.04000000000003 meter / second\n",
       "9.9    -99.24899999999981 meter  -97.02000000000004 meter / second\n",
       "10.0  -108.99999999999982 meter  -98.00000000000004 meter / second"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's position as a function of time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving figure to file figs/chap20-fig01.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3hUVf7H8ffMpEMKnUCoAQ7Se1HpltXdBUV07WDvq65lRVesiIq9IjZWXdfd/SEW7IuA0qsECBxqgNADoYZU8vvjDjhkKREycyfJ5/U882TmlplvdnE+Oeeee46nuLgYERGRcON1uwAREZGjUUCJiEhYUkCJiEhYUkCJiEhYUkCJiEhYUkCJiEhYinC7gJKMMUlAGjDCWjvOGBMFvAYMAYqAF6y1owKOvwR4CkgGpgLDrLXbQl+5iIiUpbALKGAMUD/g9WOAAVKBROBbY8xGa+0HxphWwLvAecA84BngE6B/aT/MGBMNdAU24wSgiIiEjg+ngTHXWpsXuCOsAsoYMxRIABYHbB6K0yrKBrKNMc8BNwEfAFcCX1prp/nPH+4/prm1dmUpP7Yr8HNZ/Q4iInJSegHTAjeETUAZY5oAjwCnA9/6tyXhJGt6wKHLgbb+561wWk4AWGtzjDEb/PtLG1CbAf7xj39Qt27dU/kVRETkN9qyZQtXXHEF+L+LA4VFQBljfMBHwL3W2i3GmEO7qvp/5gQcngPEBewP3Fdyf2kUAdStW5eUlJTfUraIiJSd/7nEEi6j+B4GrLX20xLb9/t/xgZsiwP2BeyP5UiB+0VEpJwKixYUcClQzxgz2P86HngD6AZswRkksdG/ryW/dvml+/cBYIyJAxpyZJegiIiUQ2ERUNbaloGvjTG/AC/5h5nvAx4xxqThdOndC7zsP/RjYJoxpi8wExgFLLTWrghZ8SIiEhTh0sV3PCOAJcBSYC4wHmcoOtbaxcC1/tdZQGvgYnfKFBGRshQWLaiSrLUdAp7nArf5H0c7djxOaImISAUSlgFVnuw/UEBufiEAXo8Hr9eDz+clwushMsKLz1ceGqkiIuFHAXUKVqzP5q+v/Uxh0bFXJfZ6PURFeImK9BET5SMmOoLYqAhiYyKIi4mgSkwkVWIjiY+LIj4ukoQq0SRUjSKxShSJVaNJqBKFx+MJ4W8lIhIeFFCnIKFKFA3qxLN7Xz7FxcUUF0PRwWKKDh6ksKiYwsIiDh4sJje/iNz8IvbsP/F7lhTh85BUNZpqCTHUTIqlRkIMNZJiqZUUS+1qcdSuHku1+Bi8XoWYiFQsCqhTULdGFV65p98x9xcXF1N0sJj8giLyCw6Sm1/ohFVeITm5hezPLSAnt4B9OQXszclnb04Be/bnsXtfPnv257FrXz77DxSQtTuXrN25rNyw66ifExnhpU71OOrWqEJyzSrUq1mFerWqUr9WVWolxSq8RKRcUkAFkcfjIcLnIcLnJS4GIPo3v0d+QRHZe/PYuTuXHXsOkLUrl6xdB9i+K4dt2QfYtjOHPfvzydy2j8xt/3t/clSkj5RaVUmpU5VGdRNoVDeeRskJ1K4Wp+ASkbCmgApzUZE+6lSPo071Y8/elJNbwNadOWzZsZ/NWfvZlLWfjdv3sXHbPrL35rFm027WbNrNr/c6Q2x0BI2TE2hSL4Gm9ZNolpJIw7oJREZoUIeIhAcFVAUQFxNJk3qJNKmX+D/79h0oIHPbXtZvcR7rtuxh3eY9ZO/NY1nGTpZl7Dx8bITPS+N6CbRokESLhtVo0bAa9WtVVUtLRFyhgKrgqsZG0rJRdVo2qn7E9l1781i7aTdrN+1m9cbdrM7cxcbt+1m1YRerNuzi6xkZAFSJjcQ0rEbLRtVo1aQGLRpVIzZa/2xEQuWxd2Yxb9nWkH1el9Pq8Mj1PU543MMPP0xubi6jR48+vO3ss89m+PDh9O9f6iX5jkvfNJVUUnw0HU1tOprah7fl5BawOnM3K9Zns2JDNivWZZO1O5cFdhsL/IsUe70emtZPpE3TGrRNrUmrpjWoGhvp1q8hIi4ZNGgQN954I3l5eURHR7No0SL27t1Lr169yuwzFFByWFxMJG2b1aRts5qHt2XtOsDydU5XYPranazZuPtwK+uzqavxeCC1fiLtm9eiXbNatGpanZgo/bMSKSulac24oXPnziQlJTFlyhTOPfdcJk6cyPnnn09kZNn9wapvEjmumkmxnJlUnzPb1wecVpZdl82SNTtYsjqLFeuzWZW5m1WZuxk/eRURPi+tmlSnk7911qRegm40FqmAPB4Pf/zjH/nqq684++yz+eabb3j99dfL9DMUUPKbxMVEHtE1mJtfSPranaSt3M6ildtZvXE3aauySFuVxbiv0qmeEE0nU4fOp9WmY4vaVFF3oEiFMWjQIAYPHsyMGTOIi4ujffv2Zfr+Cig5JTFREXQytenkD6zd+/JIW5l1+LrVzj25/Hfuev47dz0+r4c2qTXo2qou3VvXpW6NKi5XLyKnomnTpqSmpjJ69GgGDhxY5u+vgJIylVg1ml4d69OrY32Ki4vJ2LyH+cu3MW/ZVpat3cGilVksWpnFO58voXFyAt3b1KVnm2Sa1k9UV6BIOTRo0CBGjhzJq6++WubvrYCSoPF4PIfvzxrSvzl7c/KZv2wrc9K3Mm/ZVjI27yFj8x7+9cMK6lSP4/R29TijXTItGlZTWImUE8nJyXTq1ImGDRuW+XsroCRk4uOi6Nu5AX07N6CgsIjFq3Ywc8lmZi3ZzNadOUyYsooJU1ZRu1osZ7avz5kd6tEsJUlhJRKG9u7dy8aNGxk7diyXXXZZUD5DASWuiIzw0allbTq1rM3Ng9uxPGMnM9I2MT1tE9uyD/DplFV8OmUVyTWr0Ltjffp0TKFBnXi3yxYRv7Vr13L11VfTu3dvBg0aFJTP8BQXH3sto8rAGNMYWDtp0iRSUlLcLqfSO3iwmGUZO5n2y0ampW1i1968w/uapSTSt3MDenesT7X4GBerFJGykpmZyYABAwCaWGszAvepBSVhxev10LppDVo3rcH1g9qweHUWUxdsZMbiTYfvt3rvy6V0MrU5q2tDurWuQ2SEz+2yRSQIFFAStnw+Lx1a1KZDi9rcfFE75izdwuT5G1jgHxU4b9lW4uMi6dMxhXN6NDrqZLkiUn4poKRciI700atDfXp1qM+uvXlMXZjJpLnrWbtpDxOnr2Xi9LU0S0nknO6N6NMphbgY3RAsUt4poKTcSYqPZlDvVAb1TmV15i7+O2c9kxdk+rsA03jvy6X07pjCuT0a0aJhNbfLFZGTpICSci01JYnUlCSG/bE1Mxdv5rtZGSxZvYPvZ6/j+9nraNYgid+f3pgzO9TXJLYi5Yz+i5UKITrSR99OKfTtlELmtr18N2sdk+auZ9WGXbz8r19494ulnNWtIb8/o4mmWBIpJxRQUuGk1I7nuoFtuPK805j2y0a+nrGWFeud5UE+/2k1XU+ryx97NaF981q6CVgkjCmgpMKKjvQxoGtDBnRtyMoN2UyctpafFm5kTvoW5qRvoVHdeAb2TqVvpxSiIjVUXSTceN0uQCQUmjeoxt2XdeL9h8/hyt+1pFp8NOu27OXVf//CtU9+zz+/t+zel3fiNxKRkFELSiqVpPho/nS2YXC/5vz8y0Y+/2k1azbu5uPvlvN/k1YwoGtDLuibSr2aVd0uVaTSU0BJpRQZ4aV/lwb065zC4tVZTJiymnnLtvLNzAy+m5VBz3b1GNKvOc0aJLldqkilpYCSSs3j8dCuWS3aNavF+i17mDBlNVMWbGD6ok1MX7SJji1qcfFZLWjTtIYGVIiEmAJKxK9h3QTuvLQjV/yuJZ//tJrvZmWwcMV2Fq7YzmmNq3PJWS3o3LK2gkokRBRQIiXUTIrluoFtuOSsFkyctpYvf17NsoydPPbOLJo1SOLSs1rQrXVdBZVIkIVVQBlj/gA8BTQBtgHPWmvfMsZEAa8BQ4Ai4AVr7aiA8y7xn5cMTAWGWWu3hbp+qVji46K47BzDBX1S+WZGBhOmrGLVhl08+f4cmtZL5LJzDd0VVCJBEzbDzI0xycD/AX+11sYDFwMvGWM6AY8BBkgFugJDjTFX+89rBbwLDANqACuBT0L+C0iFFRsdweB+zXj7obO4YVAbqifEsGbTbka+P4e7XpzKnKVbqOzrqokEQ9gElLV2M1DLWvuNMcaLEzaFwF5gKDDSWpvtX9DqOeAm/6lXAl9aa6dZa3OB4cAZxpjmIf8lpEKLiYpgYO9U3n7wLG64oA3VE6JZs3E3T7w3m3tf+YkFdpuCSqQMhU1AAVhr9xpj4oA84HvgdWA7TtddesChy4G2/uetAvdZa3OADQH7RcpUVKSPgb1SGfvg2dwwqA1JVaNZsX4Xj4ydyfA3ppO+dofbJYpUCGF1DcovF6gCtAO+Bg74t+cEHJMDxPmfVy2xr+R+kaCIjvQxsHcq53RvxMTpaxn/40qWrtnBX1+bRpfT6nD1+adpEUWRUxB2AWWtPQjkA/OMMWOBLv5dsQGHxQH7/M/3l9hXcr9IUMVERzCkf3PO69nYPyHtKuYt28r85Vvp0ymFK393GnWq6+8lkd8qbLr4jDF9jDHzS2yOBrKBLTiDJA5pya/deumB+/xdhA05sktQJOiqxEZyxe9aMnb42Qzs3RSf18uU+Znc/PQk3v58MXv257tdoki5Ek4tqF+A+saYvwAvA92B64ALcQLqEWNMGk6X3r3+YwA+BqYZY/oCM4FRwEJr7YrQli/iSIqP5oZBbRnYK5WPvl3G1AWZfPHTGibNWc8lZ7XgD2c21ezpIqUQNi0oa+1u4HxgMLATGAtcb62dCowAlgBLgbnAeGCM/7zFwLX+11lAa5wh6iKuqlM9jnsu78xLd/elQ4ta7M8t5P2J6dzyzCSmLsjUiD+RE/BU9v9IjDGNgbWTJk0iJSXF7XKkAluwfBvvT1xKxuY9AJhG1bh+UBtaNqrucmUi7snMzGTAgAEATfy3ER0WTl18IhVap5a1ad+iFpPmrufDb5Zh12Vz3ys/07tjfYb9vjW1qpUc6yNSuYVNF59IZeDzejineyPeemAAFw9oTmSEl58WbuTmZybxz+8teQVFbpcoEjYUUCIuiIuJ5OrzWzHmrwM4s3098guK+Pi75dz6zCRmpG3S9SkRFFAirqpdPY6/Xt2VUbeeQdN6iWzLPsCov89lxFsz2bB1r9vlibhKASUSBtqk1uSFu/tw60XtiI+L5JeV27njucm89+VScnIL3C5PxBUKKJEw4fN6OO/0Jox54Cx+17MxB4uLmTBlFbc++yM//7JR3X5S6SigRMJMQpUobhvSnufv7E3zBkns2J3Lsx/OY8TYmWzK0gxeUnkooETCVPMG1Rj9597cOqQ9VWMj+WXFdm4fPZl/fm8pKNRoP6n4FFAiYczn9XBez8a8+dcB9O/SgILCg3z83XLueG4yi1dnuV2eSFApoETKgaT4aO6+rBNP3XIG9WtVZeP2/Tz4xnRe+ddC9uZoElqpmBRQIuVI22Y1efXevlx+jiHC5+WHOeu55ZlJ/LRQc/tJxaOAEilnIiN8XHZuS165py+tm9Zg9758Rn80nyfem03WrgMnfgORckIBJVJONagTz1O3nMHtF7cnLiaCuelbufXZH/lmxloOHlRrSso/BZRIOeb1eji3R2PeuL8/PdrU5UBeIW+MT+NvY2awOWu/2+WJnBIFlEgFUCMxlgeHdeOvV3chsWoUi1dncftzk/nip9VqTUm5pYASqSA8Hg9ntq/P6/f1p0/HFPILinj78yU8+OZ03eAr5ZICSqSCSawazb1Xduaha7qRFB/N0jU7uOO5KXz58xq1pqRcUUCJVFA92iTzxv396dvJaU2N/WwxD781g207c9wuTaRUFFAiFVh8XBT3XNGZB4d1JbFqFGmrnGtTP8xep/umJOwpoEQqgZ5t6/H6ff3p2TaZA3mFvPLvXxj5/hx27c1zuzSRY1JAiVQSiVWjGT60K3+5vBNVYiKYvXQLdzw3mdlLNrtdmshRKaBEKhGPx0O/zg149d7+tGtWk1378njy/Tm89p9fyM0rdLs8kSMooEQqoVrVYnniptO5bmAbIiO8fDdrHXe+MIUV67PdLk3kMAWUSCXl9Xq4oE8qL9zVh0Z149mUtZ/7X/2Z/0xaQZGGo0sYUECJVHKNkxN44a4+DOzdlKKDxXzw9TIeHjNDE8+K6yJKe6Axpg7QGagNFAFbgAXW2h1Bqk1EQiQq0scNg9rSydTmpU8Wsnh1Fnc8N5k//6kjPdsmu12eVFLHDShjTARwOXAX0B7IB7IBH1Ddf8xs4A3gE2vtwaBWKyJB1bllHV69px8v/2sh85Zt5alxczj/9MZcN7ANUZE+t8uTSuaYXXzGmD5AGnA18C7QAoiz1taz1tYBooCOwMfA7cByY0zfoFcsIkGVFB/NiOu6c8OgNkT4PHw9I4N7Xv6JDVv3ul2aVDLHa0HdA/zJWrv4aDuttcXAEv/jDWNMR+BxYEpZFykioeXxeBjYO5VWTWsw+sN5ZGzew90vTeWWwe0Y0LWh2+VJJeGp7NOdGGMaA2snTZpESkqK2+WIhJ2c3ALe/DSNKfMzAejfpQG3DG5HTHSpL2GLHFNmZiYDBgwAaGKtzQjc91sGScQBTYDokvustQtOsUYRCVNxMZH85bJOtG9Wkzc/XcyP8zawYn02DwztSqO6CW6XJxVYqYaZG2OuBLbhXJOaV+IxN2jViUhY8Hg8nNWtES/c1ZsGdeLJ3LaPe17+iUlz17tdmlRgpW1BjcIZKPECkBuMQowxZwNPA81xwnC0tfYtY0wU8BowBGd4+wvW2lEB510CPAUkA1OBYdbabcGoUaSya1Q3gRfu7M0b4xcxeX4mL32ykKVrdnDT4HZEa5SflLHSBlQC8Jq1dl0wijDGNADGA0OBz3Hut/rOGJMB9AUMkAokAt8aYzZaaz8wxrTCCc7zcFpzzwCfAP2DUaeIQEx0BHdf1ok2qTV569M0fpiznlWZu3hgaFfq1azqdnlSgZR2JokPgWFBrKMx8LG1doK19qC1di7OaMAzcEJrpLU2238B7TngJv95VwJfWmunWWtzgeHAGcaY5kGsVaTS83g8nNO9Ec/d2ZvkmlVYu2kPf3lxKjMXa2Z0KTulbUGNBhYYY64AMoAjbsi11p5Si8Va+zPw86HXxpjqQC+cYEwG0gMOXw609T9vhdNyOvQ+OcaYDf79K0+lJhE5sSb1Ennxrj68/K+FzFy8mafGzeGifs246rzT8Pk0k5qcmtIG1IfAPuArIKjrRRtjEoEvgNnAfP/mwM/MAeL8z6sepZ7A/SISZFViIxk+tCsTpqzm71+nM37yKlZl7uK+K7uQWPV/Bv2KlFppA6or0N1amxbMYowxLXCuQaUDVwCx/l2xAYfF4YQlwP4S+0ruF5EQ8Hg8DO7XjOYNknj2w3ksWpnFXS9OZfjQrrRoWM3t8qScKm0b3AJJwSzEGNMbp9X0GTDEWptrrc3GmZTWBBzakl+7/NID9/nv1WrIkV2CIhIibZvV5MW7+2AaVSNr1wEeeH0aP8wOytgqqQR+yzDzccaY14DVQEHgTmvt16dShDEmFZgIPGStfbXE7g+BR4wxaThdevcCL/v3fQxM888BONNf50Jr7YpTqUdETl7NpFhG3Xomb3++mG9mZPDKv39hZeYubhjUlsgIXZeS0ittQP3T//O5o+wrxpnd/FTcBsQDo4wxowK2vw6MAJ4HluK0+MYCYwCstYuNMdf6X9fHaYFdfIq1iMgpiozwcutF7WmeksQb49P4ZkYGGZv2MHxoV6olxLhdnpQTmotPc/GJBNWK9dmMGjeHrN251EiM4cFh3XRdSg473lx8J1pu4zcxxugGWRE5QouG1Xjh7j60alKdHbtzeeD1afw4T1MkyYkdr0P4bmPMN8aY84wxkcc6yBgTYYy5wBjzX5yFDUVEjlAtPoYnbz6D83o2pqDwIC/+cyHvfrGEooOVuwdHju+Y16CstRcYYy7EmR+vkTFmCs51oCzAA9TCWWW3J7AeeMJa+39Br1hEyqXICC+3DmlPk/qJvPVpGp9NXc36rXu578ouVI095t/AUokdd0iNf+qh9sAFwAqcMLoeuAbogjO7+UBrbXuFk4iUxnk9G/PkzaeTUCWKBcu3ce/LP7Fpu25dlP9VqlF81topaKVcESkjbVJr8sJdfXjyvdlkbN7DPS//xANXd6V9i1pulyZhRDcliIgr6lSP49k7etG9dV32HShgxNsz+XrGWrfLkjCigBIR18RGR/DgsG5c1K8ZBw8W8+b4NN6akEZR0cETnywVngJKRFzl9XoY9ofW3H1ZRyJ8XiZOW8sT780mJ7fgxCdLhaaAEpGw0L9LQ568+XTi46KYv3wb97/6M9t2BnXxBAlzpZ3qCGNMbaAdEIkzzPywU52LT0QEoHXTGjx/Z28ef3cW67bs5Z5XfuLha7tr5olKqlQtKGPMdcAG4HucNaEmBjy+DFp1IlLpJNeswug/96ZD81rs2pvH8DemMyNtk9tliQtK28V3H/A2kGit9ZZ4nOpEsSIiR6gaG8kjN/Tg7G4NyS8o4ukP5vLp5FVU9rlDK5vSBlQD4GVr7d5gFiMickiEz8sdl3Tg6vNPo7gY3p+4lDGfaoRfZVLagPoeGBDMQkRESvJ4PFw8oAX3X9mFCJ+Xr2dkMHLcHHLzCt0uTUKgtIMkFgEvGGMG4kx5lB+401p7f1kXJiJySK+O9ameGMPI92czN30rw9+czojrulMtXmtLVWSlbUH1wVkMMBZngtiuAY8uwSlNRORXrZvW4Nk7elG3RhyrNuzivld+ZqPm8KvQSjsXX79gFyIiciIpteMZfYczDH2lP6RGXNedlo2ru12aBMFvuQ+qDnA70Bqn5bUMeNtauyZItYmI/I+k+GieuuUMnvlwHvOWbeWhN6dz31Vd6NEm2e3SpIyV9j6objjXni7EWQ9qO/AHIM0Yoy4+EQmpmOgI/nZNN87t0Yj8woOMGjeHb2dmuF2WlLHStqCeB/4J3GKtPXwjgjHmNWA0oC5AEQkpn8/LbUPaUyMhho+/t7z+f4vYuSeXy84xeDyeE7+BhL3SDpLoArwYGE5+r+IMlBARCTmPx8Nl57bk9ovb4/XAP/1BpaXkK4bSBtRmoPFRtjcFdPOuiLjq3B6NeXBYN6IivHw3ax3PfDCX/IIit8uSU1TagPoQGGuMudAYk+x/DAbG+PeJiLiqe5tkHr/pdKrERjJz8WYeeXsm+w9oyY7yrLQBNRJnNol/A5nARpxrUv8BHgpOaSIiv03rpjV4+rYzqZ4Qw5LVO3jwjelk7811uyw5SaUKKGttvrX2BqAm0BPnZt0ka+291lr9iSIiYaNxcgLP3tGLejWrsGbTbv762jS2al2pcumYAWWMOd8YExnw/HzgDJyQagD0C9guIhI26lSP45nbe5GaksjmrP3c/+rPrNu8x+2y5Dc6XgtqIlAt4PmxHloPSkTCzqEbetuk1mDnnlweeH0adt1Ot8uS3+CY90FZa71Hey4iUl7ExUTy2A09efbDecxeuoW/jZnB367pTvsWtdwuTUqhtDNJ/GiMSTrK9lrGmPllX5aISNmIivQxfGhX+nVOITe/iEffmcXMxZvdLktK4ZgtKGNMX6CV/2Uf4CZjTMl7nk4DUoNTmohI2fD5vNx1aSeqxEQycfpanv5gLndd2pF+nRu4XZocx/GmOtoB3At4/I/bgMA734qBfcA9QatORKSMeL0ebrywLVViI/nXf1fw4j8XkJtXyHmnN3G7NDmG412DWowzUwTGmMnAYGttdrAL8k9MO9FaW9v/Ogp4DRiCE5AvWGtHBRx/CfAUkAxMBYZZa7cFu04RKX88Hg9XnncasdERjPsqnTfGp3Egr5DB/Zq7XZocxfG6+OKstYduHvj9oW1HOzbguJNmjPEA1wHPldj1GGBwuhITgW+NMRuttR8YY1oB7wLnAfOAZ4BPgP6nWo+IVFwX9W9ObEwEb45P4/2J6eTmF2mS2TB0vEESe40xtf3P9+HMuVfycWh7WXgMuAV4ssT2ocBIa222tTYDJ8Bu8u+7EvjSWjvNWpsLDAfOMMbozyEROa7zT2/C3Zd1OjzJ7PsT0yku1iSz4eR416D6A4duGgjFchpjrLUj/IMzAPCPHEwG0gOOWw609T9vhdNyApyWnDFmg3//yqBXLCLlWv8uDYiK9PLcR/OZMGUVefmF3HRhO7xetaTCwfGuQU092nM4fF2oHbDCWlsmt2dbazcdZXNV/8/ALsQcIC5gf8nuxcD9IiLHdWb7+kRF+nj673P5ekYGBYUHue3iDvgUUq4r7X1QzYwxU40xPfzXoeb4H+uMMT2CWN9+/8/YgG1xOF2Lh/bHcqTA/SIiJ9StVV1GXNedqEgfP8xZz8ufLKCo6KDbZVV6pZ0h4lWca00ZwFVACs7AhTeBF4JSGeAfNbjF/1mHtOTXLr/0wH3+8GzIkV2CIiIn1KFFbR69oQcxUT4mz8/k+Y8XUKiQclVpA6oXcLe1dgtwAfCVtXYl8DbQIVjF+X0IPGKMqWmMaYxzb9ahNag+BgYZY/oaY6KBUcBCa+2KINckIhVQ29SaPHZjT2KjI/j5l408++E8CgoVUm4pbUDlApHGmCo4s0p8499eF9gdjMICjACWAEuBucB4nIUSD92rda3/dRbQGrg4yPWISAXWqkkNnrz514UPn/lgrkLKJZ7SDKs0xvwDZ8n3vUAXoBHQA3gZmG6tvenYZ4c3f6ts7aRJk0hJSXG7HBEJE6s27OLht2aw70ABXVvVYfjQrkRG+Nwuq8LJzMxkwIABAE38txIdVtoW1E04w7lzgd9ba/cDXYEpwF1lVqmISJho1iCJJ28+nfi4SOamb+WpcXPJLyg68YlSZkrVggpkjEkAvNbaXcEpKbTUghKR41m7aTd/GzODPfvz6XKa05KKilRLqqyURQsKY8wt/ptgs4EdxpjNxpgHyrRSEZEw06ReIk/efDoJVaKYt2wrT42bo5ZUiJT2Pqh7gadxhpv3AnoDLwL3G2PuDF55IiLua1IvkZG3nEFClbThf1kAABJqSURBVCjmL9/GSIVUSJS2BXUbcLO19llr7Qxr7XRr7bM4c+fdHrzyRETCQ+PkhMMhtWD5Nkb9fS4FhQqpYCptQNXCGeJd0nycm3ZFRCq8wJCat2wrT/9d90kFU2kDaglHv7/oTziTt4qIVAqNkxMOj+6bk75F90kF0fFmMw80AvjKGNMTmOnf1hP4HTA4GIWJiISrJvUSeeKm0/nbmBnMXrqF0R/N4/6ruhDhK/W4MymFUv2vaa39HhgA5OHMxTcE2AN0tdZODF55IiLhKTUliSduOp0qMRHMXLyZFz/WBLNlrbQtKKy1PwE/BbEWEZFypVmDJB67sScPvzWTn37ZiM/n4c5LO2mpjjJy3CXfgZdwWkt5wATggbJa/0lEpCIwjarz6A09eGTsTCbPzyTC5+X2izto0cMycLwuvseAPwLP4iyp8Xuc2ctFRCRAqyY1GHF9j8PrSb392WItH18GjhdQQ4DLrbVPW2tH44ziG2SMiQxNaSIi5Ufb1Jo8dE03InxeJk5fy7iJ6QqpU3S8gErhyCHkc/3H1wlqRSIi5VQnU5vhQ7vi83r4dMoqPv7Oul1SuXa8gPIBh2+TttYW41yLigp2USIi5VW31nW578oueD3wyQ+W8T+udLukckuD9kVEytgZ7etx56WdABj3VTpfTV/rckXl04mGmQ8zxuwrcfyVxpiswIOstW+UeWUiIuVY/y4NyMsv5I3xaYz5NI3oSB9ndWvodlnlyvECaj3OZLCBtgDXlNhWDCigRERKOO/0JuTmF/Hel0t59d8LiY2J4Ix29dwuq9w4ZkBZaxuHsA4RkQrpwr7NOJBXyD+/tzz30Txir+1Bp5a13S6rXNA1KBGRILvsHMPA3k0pLCpm5Lg5LF2zw+2SygUFlIhIkHk8Hq4f2IazuzUkv6CIx9+dxerMXW6XFfYUUCIiIeDxeLjt4g6c0b4eObmFPPL2TDZu33fiEysxBZSISIj4vB7uubwTHVvUYve+fB5+awbbsw+4XVbYUkCJiIRQZISPB4d1o2WjamzPPsCIsTPYvS/P7bLCkgJKRCTEYqIjeOT6HjROTiBz2z4efWcWObkFbpcVdhRQIiIuqBoXxWM39qRujThWbdjFU+PmUFBYdOITKxEFlIiIS6onxPD4jaeTFB/NopVZPP+PBRQd1AzohyigRERclFyzCo/f2JMqMRFMT9vEmE/TtEyHnwJKRMRlTeol8vB1PYiK8PLtzAwt0+GngBIRCQOtm9bgvqt+XaZDM6AroEREwkaPNsncOqQDAG9NSGPaoo0uV+QuBZSISBg5t0cjrjrvNIqL4fl/LCBt1Xa3S3JNhQgoY0x7Y8xMY8x+Y8xiY0xXt2sSETlZFw9ozh/ObEJh0UFGvj+HtZt2u12SK8p9QBljooDPgX8BScBI4HtjTIKrhYmInCSPx8P1g9oenrfv0bdnsnVnjttlhVy5DyigLxBprX3JWltgrf0EWAr8yd2yREROns/r4S+XdaJtak127snjkbEz2bM/3+2yQqoiBFQrYFmJbcuBti7UIiJSZqIifTx0TTcaJyewcfs+nnh3FnkFlWe2iYoQUFWBkm3fHCDOhVpERMpUldhIHr2hBzWTYlm+LpvnPppXaWabqAgBtR+ILbEtDtBCKyJSIdRIjOXRG3pQJTaSWUu28PZniyvFbBMVIaDSAVNiW0v/dhGRCqFR3QQeuqYbET4vX01fy6eTV7ldUtBVhICaDHiMMXcbYyKNMZcC7YAJLtclIlKm2qbW5C+XdwJg3FfpTF2Q6XJFwVXuA8pamw+cB1wE7AQeAi6w1lbeu9tEpMLq1aE+1/6xNQAvfbKQxauzXK4oeCLcLqAsWGuXAGe6XYeISChc0CeVbdk5TJy2lpHvz+HZ28+kYd2Kd+tnuW9BiYhUNodu5O3ZNpn9Bwp47J1ZZO/JdbusMqeAEhEph3xeD3+5vBOmYTW2ZR/g8fdmk5tX6HZZZUoBJSJSTsVERfC3a7sfXjb+uX/Mr1D3SCmgRETKsaT4aB65vgdVYyOZvXQL736xxO2SyowCSkSknEupHc/fru1OhM/Llz+v4YufV7tdUplQQImIVACtm9bgzj85ix2++/kS5izd4nJFp04BJSJSQfTt3IDLz23JwWJ49qN5rMrc5XZJp0QBJSJSgVx6dgv6dU4hL7+IJ96dRdauA26XdNIUUCIiFYjH4+GOSzrQumkNdu7J44l3Z3OgnA4/V0CJiFQwkRE+HhzWjeSaVVizaTfPl9Ph5wooEZEKKKFKFCOu6354+Pm4iUvdLuk3U0CJiFRQKbXjGT6sKz6vh8+mrua7WRlul/SbKKBERCqwds1qcduQ9gC8OT6NtFXlZ6EHBZSISAV3dvdGXNi3GUUHixk1bi4bt5ePBccVUCIilcDQ37eie+u67DtQwOPvzGJvTr7bJZ2QAkpEpBLweT3cc0VnmtRLYFPWfp75YC6FRQfdLuu4FFAiIpVEbLQz+3lS1WgWrczinc/De2JZBZSISCVSu1ocD13TjQifl6+mr+XrGWvdLumYFFAiIpVMy8bVueMSZ2LZtyYsZtGK8BzZp4ASEamE+ndpwEX9mnHwYDFPfzCXzVn73S7pfyigREQqqavOb0XXVnXYd6CAJ96bRU5ugdslHUEBJSJSSfm8Hu69ojMN6sSzYes+Rn8UXnP2KaBERCqxuJhIHr62O/FxkcxbtpUPv053u6TDFFAiIpVccs0q/PXqrni9HsZPXsXUBZlulwQooEREBGjfvBY3DGoDwCv/WsiqDe6vxquAEhERAH5/RhPO7taQ/MKDjHx/Ntl7cl2tRwElIiKAsxrvLRe147TG1cnancuov8+loNC96ZAUUCIiclhkhI/hQ7tSMzGGZRk7GfvZYtdqUUCJiMgRqiXE8OA13YiM8PLtzAy+mZnhSh0KKBER+R/NG1Tj9oudhQ7HTkgjfe2OkNeggBIRkaPq36UhA3s3pbComFF/n8uO3QdC+vkKKBEROaZr/9Cads1qsmtvHqPGzaWgsChknx12AWWMudsY81mJbQ2NMT8YY/YaY1YbY84P2OcxxjxpjNlmjNlljHnRGBMR+spFRCoen8/L/Vd1oXa1WOz6bN4cn0ZxcWimQwqbgDLGVDXGjAaeP8ruT4A0oAZwA/CJMaapf9+NwGCgE9Ac6Ao8GPyKRUQqh8Sq0Tw4rBtREV5+mLOeb2etC8nnhk1AAV8BTYC3AjcaY1oAXYAR1tp8a+2PwBfAdf5DhgIvWWszrbXbgUeBm0JWtYhIJZCaksTt/jWkxk5IY9nanUH/zJAFlDEmyhhT9yiPOv5DLrPWDgG2lji1FbDeWhu4WMlyoG3A/vQS++oZY6oH4/cQEams+nVucHjQxNMfzGFnkGeaCGUL6nRg81EeGwGstZuOcV5VIKfEthwg7hj7Dz2PQ0REytQ1f2hNm9Qa7NyTx9NBnmkiZIMJrLVTAM9JnLofiC2xLQ7Yd4z9h4JpHyIiUqYi/IMm7n5xKssydvLeF0u4aXC7oHxWOF2DOpZ0oKExJjCEWvJrt146YErs22ytdX8qXhGRCqhafAzDh3Ylwudl4vS1zEnfEpTPCfuAstZaYBEw0hgTbYzpBwwCPvYf8iFwrzGmkTGmJs4giQ9dKVZEpJIwjapz25B2REZ4ycktDMpnlJf7hS4CxgLbgCzgOmvtEv++MUAdYAZO995/gBFuFCkiUpmc1a0RfTs3IMIXnLZO2AWUtfbRo2zbAJx3jOMPAo/4HyIiEkLBCicoB118IiJSOSmgREQkLCmgREQkLCmgREQkLCmgREQkLCmgREQkLIXdMHMX+AC2bAnOndAiInJsAd+9vpL7FFCQDHDFFVe4XYeISGWWDKwO3KCAgrlAL5yZ1UO3lrGIiIDTckrG+S4+gidUS/eKiIj8FhokISIiYUkBJSIiYUkBJSIiYUkBJSIiYUkBJSIiYUkBJSIiYUkBJSIiYUkBJSIiYUkzSZwkY0x7YAzQDlgDXGut/Z87oSsqY8zZwNNAc2AbMNpa+5a7VYWeMSYJSANGWGvHuVxOSBljkoE3gX5ALjDWWvuwu1WFjjGmB/AKYIDtwNPW2nfcrSo0jDHdgInW2tr+11HAa8AQnBl5XrDWjjrVz1EL6iT4/8/4HPgXkASMBL43xiS4WliIGGMaAOOBJ3F+/8uAUcaYc10tzB1jgPpuF+GSz3GmCKsD9ACGGmMud7ek0DDGeHF+/1estYk4/w285v/DtcIyxniMMdcD3wNRAbsewwnqVKArzr+Fq0/18xRQJ6cvEGmtfclaW2Ct/QRYCvzJ3bJCpjHwsbV2grX2oL/lOAU4w9WqQswYMxRIABa7XUuoGWO6A02BP1trc621a3H+u5jsamGhUw2oDXiMMR6gGCgE8l2tKvgeA27B+eM00FBgpLU221qbATwH3HSqH6aAOjmtgGUlti0H2rpQS8hZa3+21t586LUxpjrOhLsL3asqtIwxTYBHgGvdrsUlnXGC+VFjzEZjzGrgQmvtZpfrCglr7Q6cLq2/AwU4E50+aK0t+b1Q0Yyx1nYG5h3a4O/mTgbSA44rk+9DBdTJqQrklNiWA8S5UIurjDGJwBfAbJwujwrPGOMDPgLutdZW1oXEDv1RUoDTkhoM3FvJuvhygcuBWJzW4yPGmHPcrCvYrLWbjrK5qv9n4HdimXwfapDEydmP848yUBywz4VaXGOMaYETSunAFdbagy6XFCoPA9Za+6nbhbgoD9hjrX3U/3qRMeYdnKD62LWqQmcwcIa19j7/66nGmHdxurW+d68sV+z3/wz8TiyT70O1oE5OOs4FwUAtObKJW6EZY3rjtJo+A4ZYa3NdLimULgWGGGN2GWN24XRlvGGMecPlukJpORDnHzB0SGX6g7cBEF1iWyFOi7JSsdZmA1s48juxTL4PK9M/qLI0Gefi6N04/dAX4Qw3n+BqVSFijEkFJgIPWWtfdbueULPWtgx8bYz5BXipkg0z/wFnaPXzxph7cL6crsO5gF4ZfI8zcvVG4G2gE3ADcL2rVbnnQ5wuzjScLr97gZdP9U3VgjoJ1tp84DycYNoJPARcYK3d7mphoXMbEI/zH+i+gMczbhcmoeFvMffBuf60GfgWeNZaO97VwkLEWrsUp5vvJmAXTrfmA9baSnEd9ihGAEtwRjPPxbkNZcypvqlW1BURkbCkFpSIiIQlBZSIiIQlBZSIiIQlBZSIiIQlBZSIiIQlBZSIiIQl3agrEiTGmHE4szwfy2M4s8BPBuKttSGZKss/l+B04Gpr7YrjHOcFZgFXWWttKGoTCaQWlEjw3Ikzy3MyzmSiAN0Ctj0HzPA/33+U84Plz8Ci44UTgH9uxccpgxsuRU6GbtQVCQFjTBuc5Sma+NfLcauOGGA90N9au6SU56wGrrPWTglmbSIlqYtPxEXGmL4EdPEZY4pxVmcdjjO/3TzgSuA+4CpgDzDcWvuh//x44HmcpbaLgR+BO4+xLAI4E93uCgwnY8zDwI1ALZx1zh601n4TcM4EnNbglDL4lUVKTV18IuHnaeAunGXUGwILcIKpK/Ap8JYx5tAaPGNxguxcnLnxioHvjDHH+uPz9zjz5gFgjLnQ/1lX4sxA/RXwH2NMQsA53wJnHec9RYJCASUSfl631k621v6CM2v8PpxWjQVewFl3p4kxpilOi+hya+1cf6voKqAx8LtjvHcXnAk9D2mMs7bTOn/X4+M4k6AGLhuRjjND9RGzuIsEm/4iEgk/qwKe5wAZ1tpDF4sPrbsVDTTyP7fGHLE8WRxOq2riUd67DpAV8PojnJGGa4wx83FWR37fWnsg4Jgd/p+1f+PvIXJK1IISCT8lF7071krFEf5jOwIdAh4tgPePcc5BwHPohX+JmM44La4ZwDAgzT+o45BD3xNFpf4NRMqAAkqk/FoGRAJVrLWrrLWrcNZmGo0TUkezBWcwBADGmMHATdba7621d+K0vPYC5wecUyvgXJGQURefSDllrbXGmC+AD4wxt+GscDsSZ3DF8mOcNh9oH/DaB4w2xmzFGTHYA6jrf35IeyCbI7seRYJOLSiR8m0oTph8hrOSaSJwtrV21zGO/wpntB8A1tr/AI/gtLpWAE8Ct1trfww4pzfwrbVWXXwSUrpRV6QSMcbEARnA76y1C0pxvBdYhzNS8OcglydyBLWgRCoRa20OTmvptlKeMghYo3ASNyigRCqfF4F2psTY9JL8raeHgJtDUpVICeriExGRsKQWlIiIhCUFlIiIhCUFlIiIhCUFlIiIhCUFlIiIhKX/BwPw1Gwt0IAeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_position(results):\n",
    "    plot(results.y, label='y')\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Position (m)')\n",
    "\n",
    "plot_position(results)\n",
    "savefig('figs/chap20-fig01.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Onto the sidewalk\n",
    "\n",
    "To figure out when the penny hit the sidewalk, we can use `crossings`, which finds the times where a `Series` passes through a given value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([8.81788535])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_crossings = crossings(results.y, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example there should be just one crossing, the time when the penny hits the sidewalk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "8.81788534972054 second"
      ],
      "text/latex": [
       "$8.81788534972054\\ \\mathrm{second}$"
      ],
      "text/plain": [
       "8.81788534972054 <Unit('second')>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_sidewalk = t_crossings[0] * s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare that to the exact result.  Without air resistance, we have\n",
    "\n",
    "$v = -g t$\n",
    "\n",
    "and\n",
    "\n",
    "$y = 381 - g t^2 / 2$\n",
    "\n",
    "Setting $y=0$ and solving for $t$ yields\n",
    "\n",
    "$t = \\sqrt{\\frac{2 y_{init}}{g}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "8.817885349720552 second"
      ],
      "text/latex": [
       "$8.817885349720552\\ \\mathrm{second}$"
      ],
      "text/plain": [
       "8.817885349720552 <Unit('second')>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sqrt(2 * init.y / g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimate is accurate to about 9 decimal places."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Events\n",
    "\n",
    "Instead of running the simulation until the penny goes through the sidewalk, it would be better to detect the point where the penny hits the sidewalk and stop.  `run_ode_solver` provides exactly the tool we need, **event functions**.\n",
    "\n",
    "Here's an event function that returns the height of the penny above the sidewalk:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def event_func(state, t, system):\n",
    "    \"\"\"Return the height of the penny above the sidewalk.\n",
    "    \"\"\"\n",
    "    y, v = state\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's how we pass it to `run_ode_solver`.  The solver should run until the event function returns 0, and then terminate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>A termination event occurred.</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                             True\n",
       "message    A termination event occurred.\n",
       "dtype: object"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results, details = run_ode_solver(system, slope_func, events=event_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The message from the solver indicates the solver stopped because the event we wanted to detect happened.\n",
    "\n",
    "Here are the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>v</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>8.500000</th>\n",
       "      <td>26.97500000000022 meter</td>\n",
       "      <td>-83.29999999999998 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8.600000</th>\n",
       "      <td>18.596000000000224 meter</td>\n",
       "      <td>-84.27999999999999 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8.700000</th>\n",
       "      <td>10.119000000000225 meter</td>\n",
       "      <td>-85.25999999999999 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8.800000</th>\n",
       "      <td>1.5440000000002243 meter</td>\n",
       "      <td>-86.24 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8.817802</th>\n",
       "      <td>-4.440892098500626e-16 meter</td>\n",
       "      <td>-86.41446327683617 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                     y                                  v\n",
       "8.500000       26.97500000000022 meter  -83.29999999999998 meter / second\n",
       "8.600000      18.596000000000224 meter  -84.27999999999999 meter / second\n",
       "8.700000      10.119000000000225 meter  -85.25999999999999 meter / second\n",
       "8.800000      1.5440000000002243 meter              -86.24 meter / second\n",
       "8.817802  -4.440892098500626e-16 meter  -86.41446327683617 meter / second"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the `events` option, the solver returns the actual time steps it computed, which are not necessarily equally spaced. \n",
    "\n",
    "The last time step is when the event occurred:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "8.81780237518735 second"
      ],
      "text/latex": [
       "$8.81780237518735\\ \\mathrm{second}$"
      ],
      "text/plain": [
       "8.81780237518735 <Unit('second')>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_sidewalk = get_last_label(results) * s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is accurate to about 4 decimal places.\n",
    "\n",
    "We can also check the velocity of the penny when it hits the sidewalk:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "-86.41446327683617 meter/second"
      ],
      "text/latex": [
       "$-86.41446327683617\\ \\frac{\\mathrm{meter}}{\\mathrm{second}}$"
      ],
      "text/plain": [
       "-86.41446327683617 <Unit('meter / second')>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v_sidewalk = get_last_value(results.v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And convert to kilometers per hour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "-311.09206779661025 kilometer/hour"
      ],
      "text/latex": [
       "$-311.09206779661025\\ \\frac{\\mathrm{kilometer}}{\\mathrm{hour}}$"
      ],
      "text/plain": [
       "-311.09206779661025 <Unit('kilometer / hour')>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "km = UNITS.kilometer\n",
    "h = UNITS.hour\n",
    "v_sidewalk.to(km / h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If there were no air resistance, the penny would hit the sidewalk (or someone's head) at more than 300 km/h.\n",
    "\n",
    "So it's a good thing there is air resistance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Under the hood\n",
    "\n",
    "Here is the source code for `crossings` so you can see what's happening under the hood:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "def crossings(series, value):\n",
      "    \"\"\"Find the labels where the series passes through value.\n",
      "\n",
      "    The labels in series must be increasing numerical values.\n",
      "\n",
      "    series: Series\n",
      "    value: number\n",
      "\n",
      "    returns: sequence of labels\n",
      "    \"\"\"\n",
      "    values = magnitudes(series - value)\n",
      "    interp = InterpolatedUnivariateSpline(series.index, values)\n",
      "    return interp.roots()\n",
      "\n"
     ]
    }
   ],
   "source": [
    "source_code(crossings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [documentation of InterpolatedUnivariateSpline is here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.InterpolatedUnivariateSpline.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "**Exercise:** Here's a question from the web site [Ask an Astronomer](http://curious.astro.cornell.edu/about-us/39-our-solar-system/the-earth/other-catastrophes/57-how-long-would-it-take-the-earth-to-fall-into-the-sun-intermediate):\n",
    "\n",
    "\"If the Earth suddenly stopped orbiting the Sun, I know eventually it would be pulled in by the Sun's gravity and hit it. How long would it take the Earth to hit the Sun? I imagine it would go slowly at first and then pick up speed.\"\n",
    "\n",
    "Use `run_ode_solver` to answer this question.\n",
    "\n",
    "Here are some suggestions about how to proceed:\n",
    "\n",
    "1.  Look up the Law of Universal Gravitation and any constants you need.  I suggest you work entirely in SI units: meters, kilograms, and Newtons.\n",
    "\n",
    "2.  When the distance between the Earth and the Sun gets small, this system behaves badly, so you should use an event function to stop when the surface of Earth reaches the surface of the Sun.\n",
    "\n",
    "3. Express your answer in days, and plot the results as millions of kilometers versus days.\n",
    "\n",
    "If you read the reply by Dave Rothstein, you will see other ways to solve the problem, and a good discussion of the modeling decisions behind them.\n",
    "\n",
    "You might also be interested to know that [it's actually not that easy to get to the Sun](https://www.theatlantic.com/science/archive/2018/08/parker-solar-probe-launch-nasa/567197/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "astronomical_unit"
      ],
      "text/latex": [
       "$\\mathrm{astronomical_unit}$"
      ],
      "text/plain": [
       "<Unit('astronomical_unit')>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "N = UNITS.newton\n",
    "kg = UNITS.kilogram\n",
    "m = UNITS.meter\n",
    "AU = UNITS.astronomical_unit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>r</th>\n",
       "      <td>149597870691.0 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v</th>\n",
       "      <td>0.0 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "r    149597870691.0 meter\n",
       "v      0.0 meter / second\n",
       "dtype: object"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "r_0 = (1 * AU).to_base_units()\n",
    "v_0 = 0 * m / s\n",
    "init = State(r=r_0,\n",
    "             v=v_0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>r    149597870691.0 meter\n",
       "v      0.0 meter / s...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>G</th>\n",
       "      <td>6.674e-11 meter ** 2 * newton / kilogram ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>m1</th>\n",
       "      <td>1.989e+30 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>r_final</th>\n",
       "      <td>701879000.0 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>m2</th>\n",
       "      <td>5.972e+24 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>10000000.0 second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>dt</th>\n",
       "      <td>100000.0 second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "init       r    149597870691.0 meter\n",
       "v      0.0 meter / s...\n",
       "G              6.674e-11 meter ** 2 * newton / kilogram ** 2\n",
       "m1                                        1.989e+30 kilogram\n",
       "r_final                                    701879000.0 meter\n",
       "m2                                        5.972e+24 kilogram\n",
       "t_end                                      10000000.0 second\n",
       "dt                                           100000.0 second\n",
       "dtype: object"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "radius_earth = 6.371e6 * m\n",
    "radius_sun = 695.508e6 * m\n",
    "\n",
    "t_end = 1e7 * s\n",
    "dt = t_end / 100\n",
    "\n",
    "system = System(init=init,\n",
    "                G=6.674e-11 * N / kg**2 * m**2,\n",
    "                m1=1.989e30 * kg,\n",
    "                r_final=radius_sun + radius_earth,\n",
    "                m2=5.972e24 * kg,\n",
    "                t_end=t_end,\n",
    "                dt=dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def universal_gravitation(state, system):\n",
    "    \"\"\"Computes gravitational force.\n",
    "    \n",
    "    state: State object with distance r\n",
    "    system: System object with m1, m2, and G\n",
    "    \"\"\"\n",
    "    r, v = state\n",
    "    G, m1, m2 = system.G, system.m1, system.m2\n",
    "    \n",
    "    force = G * m1 * m2 / r**2\n",
    "    return force"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "3.5423376937972604×10<sup>22</sup> newton"
      ],
      "text/latex": [
       "$3.5423376937972604\\times 10^{22}\\ \\mathrm{newton}$"
      ],
      "text/plain": [
       "3.5423376937972604e+22 <Unit('newton')>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "universal_gravitation(init, system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def slope_func(state, t, system):\n",
    "    \"\"\"Compute derivatives of the state.\n",
    "    \n",
    "    state: position, velocity\n",
    "    t: time\n",
    "    system: System object containing `m2`\n",
    "    \n",
    "    returns: derivatives of y and v\n",
    "    \"\"\"\n",
    "    y, v = state\n",
    "    m2 = system.m2    \n",
    "\n",
    "    force = universal_gravitation(state, system)\n",
    "    dydt = v\n",
    "    dvdt = -force / m2\n",
    "    \n",
    "    return dydt, dvdt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0 <Unit('meter / second')>,\n",
       " -0.005931576848287442 <Unit('newton / kilogram')>)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "slope_func(system.init, 0, system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def event_func(state, t, system):\n",
    "    r, v = state\n",
    "    return r - system.r_final"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "148895991691.0 meter"
      ],
      "text/latex": [
       "$148895991691.0\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "148895991691.0 <Unit('meter')>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "event_func(init, 0, system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>A termination event occurred.</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                             True\n",
       "message    A termination event occurred.\n",
       "dtype: object"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "results, details = run_ode_solver(system, slope_func, events=event_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5600110.3392641265"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "t_event = get_last_label(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "64.81609188963108 day"
      ],
      "text/latex": [
       "$64.81609188963108\\ \\mathrm{day}$"
      ],
      "text/plain": [
       "64.81609188963108 <Unit('day')>"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "seconds = t_event * s\n",
    "days = seconds.to(UNITS.day)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "results.index /= 60 * 60 * 24"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "results.r /= 1e9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3iUVfbA8e+kJyQhgRBCSOhwQm+hCyrYEXt37bquq6tr+7m77mJb21rWtoqr7qprL6urqIgCovQmnRxqgAAptAQIoST5/fFOcIgkDGYmM5Ocz/PMk5n7zsx7XsEc7n3vPddVWVmJMcYYE2zCAh2AMcYYcySWoIwxxgQlS1DGGGOCkiUoY4wxQSki0AH4m4hEAwOALUB5gMMxxhhzuHCgFTBXVfd5HmjwCQonOf0Q6CCMMcbUajgwzbOhMSSoLQBvv/02aWlpgY7FGGOMh/z8fC6//HJw/6721BgSVDlAWloaGRkZgY7FGGPMkf3sFoxNkjDGGBOULEEZY4wJSpagjDHGBKWA3IMSkYHAeFVNdb+OBnYB+z3eNkNVT3Efvwh4BGcq4lTgalUtrN+ojTHG1Kd6TVAi4gKuA56sdqgnsF1VfzbNTkS6Aa8BpwPzgMeB94CR/o3WGGNMINX3EN8DwE3AX6u19wcW1vCZXwGfq+o0VS0D/ggME5HO/gvTGGNMoNX3EN84VR0rIidUa+8HpIrIYqAl8D3we1XdBHTD6TkBoKqlIrIRp9e1yt8B7913kBW526EScEGYC1wuFy4XhLlcRESEEREWRkREGOFhLiIjwoiMCCMqMtx5RIThcrn8HaYxxjQ49ZqgVHVzDYf2ANOBB4EDwHPAJ8BAIB4orfb+UiDOT2Ee5vkPFvLDwk11+o4od8KKiY4gNjqc2OgIYqIiiI2OIDYmgiYxkTSJjSQ+1vnZJDaS+JhIEuOjSIiLIrFJFFGR4T66ImOMCQ1BsVBXVe/wfC0idwBFIpKJk7xiq30kDthdH7GNzM5k776DlJdXUFkJlVRSWQkVlZWUl1dSXlHBwfJKDpZXUF5ewYHySg4eLGffgQr2HyjnwMEK9rsfu/ce+MVxxESFk9jESVZJCTEkxUeTnBjt/EyIISkxmuZNY2jeNJZoS2bGmAYgKBKUiDwIvKuqK9xNUe6fZcByQDzeGwe0cbf7XXbXlmR3bfmLP19RUcn+g+XsP1BB2f6D7N3nPMr2/fR8z96D7Ck7wJ69zmP33gPsLj3ArtL9lOzZR8me/ZTtL6ds/14Kd+wFims9Z3xs5KFklZIUS2pyLC2SY2mRHEdqchzNm8YQEW4rDIwxwS0oEhTQC8gWkcvcr58FvlDVIhF5B5jmvm81E3gU+FFVVwYm1GMTFuYiJiqCmChIbBJ19A8cQWVlJXv3HaRkz35K9uxn56597NhV5v65j5279rG9pIxtJWVsL97rJLi9B1ifv+vIMbkgJSmWtOZN3I+4Qz/TU+JpEhtZl0s2xhifCJYEdR3OfafVODF9AfwaQFWXiMi1wDigNTAbuDBAcQaEy+UiLiaSuJhI0po3qfW9FRWVlOzZz7bivWwrKWPrzr0Ubi+laMdeCneUUrhjLzt2lVG4w+mNLV699WffkZQQTesW8aSnNCEjNZ70FvFktkwgrVkc4dbzMsbUE1dlZWWgY/ArEWkHrJs0aZIVi3U7cLCCoh2l5G8rZcu2PeQfepSyeese9h848rZZEeFhZKTGk5EaT5uWCWSmJdA2LZH0FvGEh9lMRWMaq7y8PMaMGcPo0aOZMGECd911F5dcconXnx01ahRAe1XN9TwWLD0oU48iI8JIb+H0jKqrqKhkW3EZm4t2k1e0+9DPvIJdFO7YS+6WEnK3lBz2maiIMNqkJdCuVVPapSfSrlUiHVo3JSHulw1pGmN+7oFXZzFvRUG9nS+7a0vuu36w1+8vLS2lWbNmzJgxg4MHD/okBktQ5jBhYS73hIpYendpcdixvfsOkle4i40Fu9iQv4v1+btYn19C0Y69rM4rZnXe4ZM3UpvF0bF1U+eRkUTH1k1JToypz8sxxtSjMWPGEBUVRVSUb/5xagnKeC02OoLOmcl0zkw+rH333gOs31JC7uZi1m0pIXdzCeu2lFC4vZTC7aXMXPLTPmQpSbF0zkyiS5tkurRJolNGEnExNinDmKM5lt5MoKSmpvr0+7xKUCISjlOOKBtIxdlYKh9nD/maShSZRiI+NpLuHZrTvUPzQ23l5RXkFe1mTV4xazcVs2bTTtbkFbN151627tx7KGm5XJCRmkBW22S6tmtGVrtmtG4RT5jd0zIm5Pi6ak6tCUpEkoFbcernNQfWAtuAcCAFaCsiW3Bm2P1DVXf6NDoTssLDw2iblkjbtERGZmcCzv2tTUW7WblhB6s27mTlhh2s21zCxgJn2PCbORsAJ+FltWtGVrtkenRIoXNmklXSMKYRqjFBiciVOKWHJuFMA/9WVfdVe08iMBy4HFgqIveq6ht+jNeEsLAwF5ktE8hsmcCoAW0AOHCwnDWbisnJ3UFO7nZW5G5ne0kZ81YUHLohHBkRRpc2yXRr34weHVLIapdsw4LGNAK19aAGAANUtaimN6hqCc6apS9EpBVwL2AJyngtMiKcrLbNyGrbDI7vSGVlJUU795KTu53l67azbO021ueXsGztNpat3caHk1YRFuaiS2YSvTq3oFfHFLLaN7PyTsYEUEZGBqrq8++1dVAm6O0q3c8Kd7JaunYrq/OKqaj46e9tRHgYWe2S6dO5BX0llY4ZSbYuy5gQ4ZN1UCLSAegKRFc/pqr/rWOMxtQoIS6Kgd3TGNjd2c+ytOwAy9ZuY/HqrSxevZV1m4tZumYbS9ds460JOcTHRtK7cwv6Sgv6dkkltVm9FL43xviYt7P47gYew9kVqXpJ7krqaesLYwDiYiIZ0C2NAd2chLWrdD9LVm9l4coiFmghBdtLmb54M9MXO7u7ZLaMp39WS7KzWtKtQ3MiI6xckzGhwNse1N0495f+pqoVfozHmGOWEBfF0F7pDO2VDsCWrXtYuLKQH1cWsWhVERsLdrOxYDefTl1DbHQ4vTu3ILtrSwZ0S6OZLRw2Jmh5m6BcwH8tOZlQ0CqlCa1S2nP60PYcLK8gJ3c781YUMD+nkNwtJcxams+spfnAIrq0SWJg9zQGdW9F27QE2/3YmCDibYL6B/AHEfmNqu73Z0DG+FJEeBg9OqbQo2MKV5/ZnaIde5mfU8Dc5QUsXFXEyg07WblhJ299lUNqszgGdU9jSM9WdGvf3CZaGBNg3iaoD4HvgctEJB84rCelqh18HZgx/tAiOZbThrTjtCHtKNt/kEUri5i9LJ+5ywso3F7K5z+s5fMf1tI0PorBPVoxpGcrenVqYfetjAkAbxPUW0Au8B5Q6rdojKlHMVERDOrRikE9WlFRUcnKDTuYtXQLMxZvYcu2PXw9az1fz1pPk5gIBnZP47g+renbJdWSlTH1xNsElQX0DpVdbI05VmFhLnd5pWZcNbobuVtKmLlkCzMWb2Z9/i6mzM9jyvw8msRGMqRHK4b3aU2vzilE2AaOxviNtwlqLiCAJSjT4LlcLtqnN6V9elMuOzWLTUW7mbZoE9MWbiZ3Swnfzt3At3M3kBAXxbDe6ZzQL4Ou7ZpZgVtjfMzbBPU28C8ReRdYQ7W1UKr6oq8DMyZYtG4Rz8UnCRefJGws2MUPCzfxw8JN5BXuZsLMXCbMzCU1OZbj+2VwQr8M2qQlBjpkYxoEbxPUH4DdwJgjHKsELEGZRiGzZQKXnZrFpacIuVtKmLogj6kL8ijcsZcPJ63iw0mr6JDelBOzMzmhXwZJCT8rvGKM8ZJXCUpV2/s7EGNCiecw4JVndGPZum18Nz+P6Ys2sXZzMWs/K+b18cvI7tqSUQMyye6aZpMrjDlG3pY6ulBVPzxCe3fgZVU9zueRGRMiwsJc9OyYQs+OKdx4bk/mrihg8tyNzMspYPayfGYvyychLooT+mdw8sA2tE9vGuiQjQkJ3g7xvSki0ar6FoCIRAP3AXcBs/wVnDGhJioynGG90hnWK50du8qYuiCPb+dsYH3+rkNrrDplJnHqoLaM6Nva9rUyphbeJqiLgXdEJBZYD7wENAFusA0KjTmy5IQYzjm+E2eP6MiaTcV8O2cD3y3IY/XGnazeuJPXPlvK8D6tOXVwW7q0SbYyS8ZU4+09qM9E5Azgf0ACToL6s6oW+zM4YxoCl8tFp4wkOmUkcc2Y7kxftJmJs9ezbO02vpmzgW/mbKBDelPOGNaO4/tmEBPt9S44xjRotW353q1a01bgd8ArQCHQWkRaA6jqcr9FaEwDEh0ZzsjsTEZmZ5JXuIuJszcwae4G1m4u5oUPF/Gvz5cxMjuTM4a2J7NlQqDDNSagavun2lKcKeRV4w5VW5i6gAfcj6r2Y9pvW0QGAuNVNdX9OhV4Fhjl/v6vgNtUdYf7+JvARcBBj6/ppaprj+W8xgSTjNQErh3TnV+dlsX0xZv5akYuK3K3M37aOsZPW0fPjimceVx7BvVoZYVrTaNUW4Ly+dRyEXEB1wFPVjv0KlDsPmck8B+cCuqXuY/3A85R1Qm+jsmYQIuKDOfE/pmc2D+TdZuL+XJGLt/N38iSNVtZsmYrqcmxjB7WgVMGtSE+LirQ4RpTb2pMUKq63g/newAYDfwV+DOAiIThVEd/QFX3uNteAV5wP4/FqQW40A/xGBNU2qc35eYLenP16G5MmreB8dPWsWXrHv49fhnvTMxhZP9MxgzvYMN/plGo77ux41R1rIicUNXg3gTxnGrvOwf40f28D87Q3isiMhjYCIxV1fH1EK8xAdEkNpKzhnfkzGEdmJdTwOc/rGXhyiK+mpnLVzNzye7aknNP6EjPjik2+880WPWaoFR189HeIyJ34SSooe6mBOAHnN7XIuAs4AMRGaKqi/wVqzHBICzMxcBuaQzslsaG/BI+n7aOyfM2Mm9FAfNWFNChdVPOPaETx/VOt8rqpsEJmvmsIhIJPI9T72+kquYAqOpEYKLHWz8WkWtwEpUlKNNotElL5OYLevOr07L4amYu46etZe2mYp56ez5vfLGcs0d04NTB7Yi1aeqmgQiKv8kikgB8jtNbGqiqmzyOjQGaVVsQHAWU1W+UxgSHpvHRXHKycN4JnZgyP49Pp64mr3A3r322jPe/WcmZx3VgzPAOJDaxCRUmtHlbi68zzsy7/jiz7A4b9K6aLl4H7wFhwHBVrb5jbzjwrIisAObjVLUYClxfx3MaE9KiIsM5dXBbTh7YhnkrCvho8ipW5G7nvW+UT6au5tTBbTn3+E6kJMUGOlRjfhFve1AvA62Ax4ESXwYgIr2AM4B9QKGIVB3aqaoZqvqpiNwLvAukATnAmaq6wZdxGBOqwsJcDOyexoBuLVm2dhsfTl7FgpxCPvt+LV9OX8fI7DZcOKozac2bBDpUY46Jq7Ky8qhvEpFdwAmqOt//IfmWiLQD1k2aNImMjIxAh2NMvViTt5OPJq9ixuLNVFQ6SWxk/0wuPKkz6SnxgQ7PmEPy8vIYNWoUQHtVzfU85m0PaiMQ4+O4jDF+0jEjiXuuHMDmot28/+1KvluQx7dzNzB5/kZO6JfBRSd1oXULS1QmuHmboP4CvCQiDwGrgP2eB60WnzHBKb1FPLdf2o9LThY+nLSSyfM2MnneRr6bv5ER/TK49BSxHpUJWt4mqKrNCt8/wrFjrsVnjKlfrVKacOvFfbnopC58NHmVs/XH/Dy+/3ETJw9sw0UndSE1OS7QYRpzGG8TlG35bkwDkNa8Cbdc2IcLR3XhvYnK5Hkb+HrWeibN3chpQ9py0aguJCfaaL4JDt7uB7UeQEQ6Ad1wpoSvUFX1Y2zGGD9p2SyO2y7py/kjO/Hu18r3Czcxfto6Js7ewFnDO3D+yM7Ex9puvyawvF0HFQ+8BlyIc//JBUSIyDfA+VVFXo0xoSUjNYG7r8jmglGdeXtCDrOX5fPR5FV8PSuXi07qwhlD2xMVaSP4JjC8Ld71NNATGAzE4szoGwKk46yNMsaEsPbpTfnztYN48tbh9OjYnF2lB3jts2X85vFJTJ63gfKKoy9HMcbXvL0HdR5wrqrO8WibIyI3Ax8Bt/g8MmNMvZO2zXjkpmHMzynkjS+Wk7ulhL+/+yOffLeGa87sTr+suhaNMcZ73iaoMJwt36vbDtgcVWMaEJfLRXbXlvSVVKYu2MhbE3LI3VLCfa/MpF9WKteN6U6btMRAh2kaAW+H+L4H7heRQ9UnRSQauA9nKwxjTAMTHuZiZHYbxt0zimvO7EZcTAQLcgr53VPf8eLHiyjevS/QIZoGztse1F3ANGCjiFTtbNsbp6L4af4IzBgTHKIiwznvxM6MGtCGt7/O4euZuXw1I5epC/K4aFQXzhrRgcgIm0hhfM+rHpSqrga6Ag8Ba4ClOFu2d63at8kY07A1jY/mt+f35rm7TqR/ViqlZQd5/Yvl3PLEFOatKAh0eKYB8no/KFXdAbzgx1iMMSGgbVoi998whAU5hbz62RI2FuzmgVdnkd21JTec08NKJxmfqTFBiUgh0E1Vt4pIEU5JoyPywX5QxpgQ0y8rlec6n8j4aet4d2IO81YUsHBlEecc35GLTupiO/uaOqvtb9DdwC6P57YQwhhzmIjwMM45viPH923NG18uZ9LcjXw0eRVT5m/khrN7MrRXK1wu19G/yJgj8Go/qFBm+0EZU390/XbGfbKE1Rt3Ak4v6zfn9qJVim2WaI7sF+0HJSJ/8/YEqvp/vzg6Y0yDIW2b8eStI5g4K5c3vlzBgpxCbnliMhed1IXzTuxks/3MMaltiG+Al9/RsLtgxphjEh7m4vSh7RncsxX/+nwZ383P460JOUyZn8dN5/eid+cWgQ7RhIgaE5SqnlifgRhjGpbkhBjuvKw/Jw9sw0sfLyavcDd/HjeDkwe24dox3YmPizr6l5hGrbYhvjO8/I5KVf3KR/EYYxqYXp1a8NydJ/Lf71bx/jcr+WbOBuauKODGc3syrFe6TaIwNaptiG+8l99hO+oaY2oVGRHGxScJw3ql88KHi1i2dhuPvzmPQd3TuOn8XjRvGhvoEE0Qqm2Iz9s6fcYY45WM1AQeuWkYX89ez+vjlzF7WT5L1mzl6jO7c9rgttabMoepbYivG5CjqhXu5zWpVNUVvg/NGNMQhYW5OH1IOwZ2a8lLHy9m9rJ8XvxoETMWbeZ3F/UhtVlcoEM0QaK2XtJSIMXj+RL3zyM9jDHmmDRvGsu91wzk/67IJiEuioWrirjlySl8PSuXhr4+03intntQ7YEij+fGGONTLpeL4X1a06Njc176eDEzl2zhhQ8XMWPxFm65sA8tku3eVGMWkEoSIjIQGF9Vw8+9z9QLwAVAOfC0qj7q8f6LgEeAVsBU4GpVLfTyXO2wShLGBL3Kykp+WLiJcf9dzK7SA8TFRHDD2T0YNaCN3ZtqwH5RJQlPItIHeAboDkRXP66qXm2vKSIu4DrgyWqHHgAE6Ag0BSaIyCZVfdN9/+s14HRgHvA48B4w0ptzGmNCg8vlYkTfDHp2TOEfHy1i9rJ8nn1/IXOWF3DzBb1pGv+zXz2mgfO23PAbwE6cjQvL6nC+B4DRwF9x9pOqchVOr2gHsENEngRuBN4EfgV8rqrTAETkj+73dFbVVXWIxRgThJITY7j3moFMmb+Rcf9dwswlW8jJ3c7vL+lHvyzbOKEx8TZBdQKyfTBbb5yqjhWRE6oaRCQJZ+huucf7coCe7ufdcHpOAKhqqYhsdB+3BGVMA+RyOdvNd++QwtPvzGf5uu3c98pMzhzWnqvHdCc60pZeNgbernX6AehV15Op6uYjNFftblbq0VYKxHkcL+VwnseNMQ1Uy2ZxPPLb47jyjK6Eh7kYP30dt//9O9ZuKg50aKYeeNuD+jUwS0RGA2uBCs+DqvpgHWLY4/7pOV0nDtjtcbz6VB7P48aYBiw8zMWFo7rQV1J56u35bCzYzZ3Pfs+1Y7pz5nHtbQJFA+ZtD+p+IBXoj3MPaYzH48y6BOC+75SPM0miShY/Dfkt9zwmInFAGw4fEjTGNHCdMpJ45o4TOH1IOw6WV/DPT5fw8L/nULJnf6BDM37ibQ/qIuAcVfW2Pt+x+g9wn4gsxhnSuwt41n3sHWCa+77VTOBR4EdVXemnWIwxQSo6MpzfXtCb3l1a8Pz7PzJ7WT63PTWFOy/vT4+OKUf/AhNSvO1BbccZ2vOXsTgVKZYBc4GPgXEAqroEuNb9eivOVPcL/RiLMSbIDeuVzrN3nkhW22S2Fpdx70vTeXeiUl5hFSgaEq8W6orIxTj3oe4G1gAHPI+ravVJDEHDFuoa03AdLK/gna9z+GjyKioroXfnFO66PJukBFszFSpqW6jrbQ/qBWAETu9mO7Cr2sMYY+pdRHgYV57RjQd/PYSk+GgWrdrKbU9/x/J12wIdmvEBb+9BXeDXKIwxpg76dEnlmTuO52//mcfyddv504vTufrM7pw9ooPN8gthtW23kaSqOwFUdao3XyYiye5ZecYYU6+aN43l4ZuG8cYXy/l06hpe+2wpy9dt47aL+9IkNjLQ4ZlfoLYhvqkico+IND3al4hIioj8Gfjed6EZY8yxiQgP47qzevDHqwYQFxPBzCVbuP2ZqazfUhLo0MwvUNsQ3zCcmnl5IjIdmIAzy24r4AJaAL2B44HhOPX6hvk1WmOM8cLQXum0a5XIo2/MJXdLCXc99z23X9qPob3SAx2aOQY19qBUdbeq/h7ojLP+6DLgS2A+Tm28z4DzgFlAlqr+TlXtnynGmKCQ3iKeJ24dzoi+rSnbX86jb8zlra9WUGFT0UPGUSdJqGo+ThXyB0QkDGgOVKiqTZMxxgS1mKgI7rq8P50yknh9/DLe/3YlazcXc+dl/e2+VAjwdpo5AKpaoapFlpyMMaHC5XJx7gmduO+GIcTHRjJ3eQF3Pvs9eYW2QibYHVOCMsaYUNVPUnn698fTrlUim4qcgrPzcwoCHZaphSUoY0yj0SqlCX/73XCG9mpFadlBHnx1Fl9M82cVN1MXlqCMMY1KbHQE91wxgItP6kJFJYz7ZAkv/3cx5eUVR/+wqVfeVpIAQEQi3Z85bGl2MNfiM8aY6sLCXPzq9K60To3nufcXMn76OjZv3cP/XZFtkyeCiFc9KBEZLCKLgDKcjQKtFp8xJuSd2D+Th28aSmKTKBZoIXc//wP52/Yc/YOmXng7xPcMUAycA4w8wsMYY0JSt/bNeeq2EWS2TGBjwS7ueu57Vm6wim3BwNshvp7AYPfeTMYY06CkNW/CE78bzmNvzmXhyiL+9NJ07rkimwHd0gIdWqPmbQ9qBWA1QowxDVaT2Ejuu34wI7Mz2be/nL/+ew5fz1of6LAaNW97UM8Dr4jI88AqYL/nQVX90teBGWNMfYsID+P3l/QlJSmWD75dyQsfLmRb8V4uPUVs244A8DZB/dv98/EjHKsEwn0TjjHGBJbL5eKK07uSkhTLuI8X8e5EZevOvfz2gt5EhNvKnPrkVYJSVftTMcY0KqcPaUezhGj+9tZ8vpmzgR279nHPldnERB3T6hxTB5Z4jDGmBoN6tOLhm4aSEBfFvBUF3P/KLPbsPRDosBoNr/4pICJFOEN5R6SqqT6LyBhjgkhW22Y8dvMwxv5zJsvWbuNPL03nwV8PoWl8dKBDa/C87UHdBdzt8fgj8Cpw0P3aGGMarDZpiTx+y3BaNW/C2k3F3PPCNIp27A10WA2et/eg3jhSu4jMBW7D2U3XGGMarJbN4nj8luMY+8+Z5G4p4Z5//MBDNw6ldYv4QIfWYNX1HtQiYJAvAjHGmGCXnBjDo78dRlbbZIp27OUPL0xj3ebiQIfVYHl7D6rbEZqbAn8GVtc1CBG5HHi5WnMsMAkYg1Pvz3Pt1QxVPaWu5zXGmGMVHxfFgzcO5ZF/z2HhqiL+9OJ0HrpxKJ0ykwIdWoPj7XzJpTiTJKqvVNsIXF3XIFT1beDtqtci0heYiHN/qyewXVWt5ogxJijERkcw9vpBPP7mPGYvy+fPL8/gwV8PoUub5ECH1qB4O8TXHujg/tkeaIdT+qidqk7xZUDuLT3eBu5X1UVAf2ChL89hjDF1FRkRzj1XDmBIz1bs2XuAv7w8A12/PdBhNSheJShVXQ9sA/Ldz+OBXwEj/BDTzcBe4EX3635AqogsFpECEflQRFr74bzGGHNMIiPC+L8rsg/t0PuXl2eSk2tJyle83Q/qNGAzMNydHH4AbgEmiMi1vgpGRKJwhvXuV9WqdVd7gOnAKEBwktcnvjqnMcbURUR4GHf/Kpvjeqezd99Bxv5zBsvXbQt0WA2Ct0N8jwDPAlOBa4DtQEec+0/3+DCe04AK4IuqBlW9Q1V/p6pFqroTuAMYICKZPjyvMcb8YhHhYdx1eX9G9G3N3n3l3Ode1GvqxtsE1RV4RVUPAGcB41W1HJgJtPFhPGcDH6hqRVWDiDwoIl093hPl/lnmw/MaY0ydhIeHccel/TihXwZl+8t54NWZtvFhHXmboIqANiLSFmfSQlUPJxtn6M9XBuMM53nqBTwlIkkikoTTk/tCVYt8eF5jjKmz8PAwfn9pP0b0cXpSY/8509ZJ1YG3CeoV4FOc5LEImCQiNwFvAv/wYTzt+HnCuw7YgbPeKhdnPdQVPjynMcb4THiYi9sv68eg7mmHZvdtLNgV6LBCkquyssYasIcRkXOBtsBbqrrV/TpcVT/yZ4B1JSLtgHWTJk0iIyMj0OEYYxqJAwfLeei12fy4sohmiTE8fstxpDVvEuiwgk5eXh6jRo0CaK+quZ7HvN7YRFU/qe21McaYn0RGhPOnawZy/yuzWLZ2G/eOm8Fjvz2OFsmxgQ4tZNh+UMYY4ycxURGMvW4Q0iaZwu2l/OXl6ezYZfO7vGUJyhhj/CguJpL7bxhMh/SmbCraw/2vzKK0zDY99IYlKGOM8bP4uCge+PUQWqU4+0k98vocDhysOPoHGzlLUMYYUw+SEqJ58NdDSEqIZtGqrTzz3gIqKrybpNZYebvdRh/gGaA78LN9jlU10cdxGWNMg5PWvAn3Xz+YP744ne9/3ERyQuM9EbIAABnCSURBVAzXndUdl6v6RhEGvJ/F9wawE2frd7vDZ4wxv1DHjCT+dPUAHnh1Fv/7fg3NEmM478ROgQ4rKHmboDoB2aq6wp/BGGNMY9CnSyq3X9qPJ96az7/HLyMpIZqR2VZetDpv70H9gFNyyBhjjA+M6JvB9Wf3AOC5939k4crCAEcUfLztQf0amCUio4G1OBXHD1HVB30dmDHGNHRnj+jItuIyPvluNY+9MZcnbh1BZsuEQIcVNLztQd0PpOIUih0NjPF4nOmXyIwxphG4enQ3BvdIY0/ZQR56bTYle/YHOqSg4W0P6iLgHFUd789gjDGmsQkLc3HnZf255x/TWLupmEffmMODvx5KZIStAvL2v8B2nKE9Y4wxPhYTHcFfrh1Es8Rolq7ZxksfL8LbQt4Nmbc9qLuB50XkbmANcFidDlUt9XVgxhjTmKQkxfLnawfxh39M55s5G8hIjee8EzsHOqyA8rYH9QIwApiL05vaVe1hjDGmjjpnJnPHpf0AeP2L5cxauiXAEQWWtwnqAuAkYGQND2OMMT4wrHc6V5zelcpKePLt+eRuKQl0SAHj1RCfqk4FEJFYoDNOYlujqtZ7MsYYH7twVGc2Fu7iu/l5PPL6HJ7+/fHEx0YGOqx651UPSkTCReRxnK3XfwQWAEUiMk5EvN700BhjzNG5XC5uvqA37dMT2bJ1D39/p3EWlvV2iO9h4FfAlUCm+3ElcAYw1j+hGWNM4xUTFcGfrh5IfGwkc5bn8+HklYEOqd55m6CuBG5Q1Q9UdbOqblLVD4AbgWv8F54xxjReac2bcOfl/XG54O0JOSzIaVzlkLxNUPHA6iO0rwVSfBeOMcYYT9ldW3LpKVnuSRPzKNjeeFb1eJug5gI3H6H9FmC+78IxxhhT3cUndSG7a0t2lR7gkdfnsO9AeaBDqhfeTnC4B/hORE4AZrnbBgPtgNN8H5YxxpgqTjmkftz+zFTWbirmpY8X8ftL+gU6LL/zqgelqvOAvsBEnAkSKcDnQJaqzvRfeMYYYwDi46L409UDiYoMZ9LcjUyZvzHQIfmdt1u+jwfuVtW7/RyPMcaYGrRPb8qN5/bk+Q8W8tLHi8hq24xWKU0CHZbfeDvENxjwaw14EbkWeBnY59F8M/AuTqmlC4By4GlVfdSfsRhjTLA6eWAbFmgh0xdt5sm35/H4LcOJCG+Ylc+9TVB/B/4jIn/Hmbm31/Ogqi73QSz9gKdU9Q+ejSLyKCBAR6ApMEFENqnqmz44pzHGhBSXy8UtF/Zh5YYdrNywk7cn5HDV6G6BDssvvE27D+H0ot7HmdG3FFji8dMX+gMLj9B+FfCwqu5Q1VzgSZz1V8YY0yjFx0Zy1+X9CXPBx1NWsWhVUaBD8osaE5SInC4i0e6X7Y/w6ODxs05EJBzoBVwhIptFZLWI/EFEkoFWgGcPLQfoWddzGmNMKOvWvjmXnCxUVsLT78ynePe+o38oxNQ2xPchkAXkAVOAAaq6zU9xtADmAW8A5wFdgf8BUe7jnivTSoE4P8VhjDEh46KTurBwVRHL123nufcX8udrB+JyuQIdls/UlqB2AmNFZDrOeqdLReSIdd/rej9IVfOB4z2aForI88Dp7texHsfigN11OZ8xxjQE4eFh3Hl5f2596jvmLM/ny+nrGH1cnQe1gkZtCer3wCPAOUAlcL/7Z3WVQJ0SlIh0By5S1fs8mqOAMiAfZ5LEJnd7FocP+RljTKOVmhzH7y7sw2NvzuVfny+jr6SS3iI+0GH5RI0JSlU/Aj4CEJEKoJuq+qtS4U7gThHJA17DWRR8K04ppWXAfSKyGKcm4F3As36KwxhjQs6w3umc2D+DKfPzeO6DhTxy0zDCwkJ/qM/bShJhfkxOqOom4Cyc2XklwMfAQ+4kORZntuAynBmEHwPj/BWLMcaEohvO6UlSQjTL1m7jyxnrAh2OTwTNZoOqOhnIPkJ7Gc6C3SMVqzXGGAMkxEXx2/N78cjrc3nji+Vkd21JWvPQrjLRMJcfG2NMIzSkZzrH9U6nbH85L3y4kMrK0N6F1xKUMcY0IL85rxeJTaJYtGorE2evD3Q4dXJMCUpE4kWkr4hEi0iCv4IyxhjzyzSNj+bGc51aBq99toyiHXuP8ong5VWCEpEoEXkRZ7bdXKA18C8RGS8iTf0ZoDHGmGMzvE9rBvdIY+++g7zwUegO9R1LLb6hwHCctUkAT+As4H3a92EZY4z5pVwuFzed35v42EgW5BQyaW5o7h3lbYK6CLjFvTlhJYCqzgFuAMb4KTZjjDG/ULPEGG44pwcA//p8GbtK/bpjkl94m6BScSo6VFeC1cUzxpigdGL/THp2TGFX6X7emZAT6HCOmbcJ6nvgNo/XlSISBfwFmObzqIwxxtSZy+Xi1+f2JCzMxZcz1rFuc3GgQzom3iaoW4HRIpIDxACvA7nAccDtfonMGGNMnbVrlcgZQ9tRUQkvf7IkpCZMeFvqaBXOFhiPAc8AC3CKx3ZR1RV+i84YY0ydXX5qFolNoli2dhs/LNx09A8EiWNZB3UcsFFV71bV24HOwAD/hGWMMcZX4uOiuPKMroAzYaJs38EAR+Qdb9dBXQ+Mx0lKVZKAr0TkUn8EZowxxndOGtiWThlN2VZcxgeTVgY6HK9424O6B7hGVQ9VEVfVG4DrcKqNG2OMCWLhYS5uPK8XAJ98t4bNW4N/31dvE1Q6zpbs1c3BWaxrjDEmyGW1bcbI7EwOllfw2v+WBTqco/I2QS3C6S1VdyW2u60xxoSMq0d3IzY6gjnL85m3oiDQ4dTK2/2g7sW533QyTk+qEugH9MQqSRhjTMhITozh0lOEf32+jDe+WE4/SQ3a3Xe9nWY+BegDTAEygJbu51mqOsl/4RljjPG1M49rT0rTGHK3lDB98eZAh1Mjr3fUVdUc4G4/xmKMMaYeREaEc9HJwosfLeLdiTkM7ZVOeBD2orxKUCLSAvgz0B+IBA67ElUd6PvQjDHG+MtJA9rw0eRVbCzYzfc/5nFi/8xAh/Qz3k6SeA24BJgNfAl8Ue1hjDEmhERGhHHpyV0AeHeiUl5eEeCIfs7bIb4RwNmqOtWfwRhjjKk/J/bP5MNJq9i8dQ+T523k5EFtAx3SYbztQe0AQqsMrjHGmFqFh4dx6SkCwHvfruTAweDqRXnbgxoLvCAitwOrgMN2vlLVUl8HZowxxv+G983gg0kr2Viwm2/nbuD0Ie0CHdIh3vagngYGAbOAbcCuag9jjDEhKDzMxaWnZAHwwTfK/gPlAY7oJ972oC7waxTGGGMCZlivdNq1SiR3Swlfz1rPmOEdAh0S4GWCqm1yhHtn3TpzV6l4DKdieiHwhKq+LCLROL00z2HFGap6ii/Oa4wxjV1YmIvLTs3ikdfn8OGklZw8qA0xUV4vk/Ubb9dBpeKUO+oGhLubXUC0uy2pLkGISCbwMXAV8D+c9VZfi0guzpDidlVNq8s5jDHG1GxwjzQ6ZTRldV4xE2et56wRHQMdktf3oF4BzgFWAMOBJUApMBh4wAdxtAPeUdVPVLVCVecC3wHDcJLVQh+cwxhjTA1cLhcXneSsi/pi+joqKgK/Nby3CeoE4CpVvRVYCvxHVUcDj+KskaoTVf1BVX9T9VpEmuEkwh9xitKmishiESkQkQ9FpHVdz2mMMeZwA7ulkZIUy+ate1i4qijQ4XidoKKBNe7nK3CSBsDrwFBfBiQiTYHPcKpW/A/YA0wHRgEC7AU+8eU5jTHGOOuiThviLNb9cvq6AEfjfYJSnB4NOAlqiPt5PBDrq2BEpAvOVPYC4AL3cN8dqvo7VS1S1Z3AHcAA930rY4wxPnTKoLZEhLuYuzyfwh2BXeLqbYL6G/BvEbkCeB+4VEReB94GvvdFICIyAqfX9ClOcipztz8oIl093lo1a7DMF+c1xhjzk+SEGIb1ak1FJUyYmRvQWLzdD+ptnPtQS1R1JXAmEAdMBa6paxAi0hEYD4xV1T+qqufduV7AUyKSJCJJwLPAF6oa+AFSY4xpgM4Y1g6AibPXc+Bg4BbuepWgRGQssEhVFwKo6reqehFwD/BHH8RxM5AAPCoiuz0ej+NsNb8DWA3k4qyHusIH5zTGGHMEXds1o316IsW79zNtUeA2NKxxHZR7plxT98v7gMkisr3a2/oAv8G5L/SLqeodR/mOy+vy/cYYY7zncrkYPaw9L3y4iC+mrwvYXlG19aAG4EwpX+J+/b37tefjLZz7UMYYYxqQ4/tm0CQmAl2/g9V5OwMSQ40JSlU/xVlA2xGnasRAoL3Hox2Qoqo3+D1KY4wx9SomOoJRA9sAgZtyXmupI1Xd4H56WCJz19/rBQQmrRpjjPG7M4a257Pv1zL1x01cO6Y78XE+Kb3qNW8nSXQUkakiMlhE4oA57sd6ERns1wiNMcYEROsW8fTp0oL9B8r5du6Go3/Ax7xdB/UCTkXxXJwZdBk4VR1ewtkryhhjTAM0elh7AL6cnlvv9fm8TVDDgdtVNR+naOwXqroKp4hsH38FZ4wxJrAGdEujRXIsW7btYfHq+l1+6m2CKgMiRaQJcDzwlbs9DSj2R2DGGGMCLzzMdWia+ZzlBfV6bm8T1Nc4vaWPcbbZ+FxERrnbPvNTbMYYY4JA/6xUAOatCM4EdSMwD6cnNVpV9+Csk/oOuN0/oRljjAkG0iaZ+NhItmzdw+ai3fV2Xm+3fN8N3Fat7TG/RGSMMSaohIeH0VdS+WHhJublFHBWi/h6OW9tpY4+AK5X1RL38xq56/IZY4xpoLK7Oglqfk4hZw2vn+3ga+tB7QEqPZ4bY4xppPqKcx9q6eqtlO0/SEyUVwNwdVLjGVT1miM9N8YY0/gkJ8TQKaMpq/OKWbpmG9ldW/r9nEdNgSKSgrP/U3cgEWda+ULgS/cOt8YYYxqB/l1bsjqvmPkrCuolQdU6i09E7gDW41SSOA3oBpwK/AvYKCK3+j1CY4wxQSE7y0lK83IKqKz0f1WJ2iZJXAM8BNwFvK6qez2OxQBXA0+IyCZV/djfgRpjjAmszm2SSYiLJH9bKZu37qG1n2fz1daDuhW4R1Vf8kxOAKpapqrjgAeoNv3cGGNMwxQe5jo0WWJ+PSzarS1BdcGpIFGbz4CuvgvHGGNMMOtfNcwX4AQVC5Qc5fPFQDPfhWOMMSaY9ZNUXC5YunYbZfsO+vVcRyt1VL+11Y0xxgS1pIRoOmUkceBgBYvXbPXruY42zfxqEamt8FKCL4MxxhgT/PpntWTVxp3MX1HAwG5pfjtPbQlqA3CTF99R/9ssGmOMCZjsrqm8940yL6eQyspKXC6XX85TWyWJdn45ozHGmJDWKTOZhLgoCreXkle4m8yW/hlM83a7DWOMMQZwpptX7RE1P6fQb+cJiQQlIr1FZKaI7BGRJSIyINAxGWNMY3YoQflxunnQJygRiQL+B7wPJAEPAxNFJDGggRljTCPW12O6+V4/TTcP+gQFnABEquozqnpAVd8DlgEXBzYsY4xpvJrGR5PZMoGD5RVs8tMuu/7f0KPuugErqrXlAD0DEIsxxhi3oT3TqajY5LeafKGQoOKB0mptpUBcAGIxxhjjdvlpWVx+Wpbfvj8Uhvj24JRd8hQH+KdPaYwxJiiEQoJaDki1tix3uzHGmAYqFIb4pgAuEbkdZ+PE84FewCcBjcoYY4xfBX0PSlX3A6fjJKbtwL3AOapaFNDAjDHG+FUo9KBQ1aXAcYGOwxhjTP0J+h6UMcaYxskSlDHGmKAUEkN8dRQOkJ+fH+g4jDHGVOPxuzm8+rHGkKBaAVx++eWBjsMYY0zNWgFrPBsaQ4KaCwwHtgDlAY7FGGPM4cJxktPc6gdclZWV9R+OMcYYcxQ2ScIYY0xQsgRljDEmKFmCMsYYE5QsQRljjAlKlqCMMcYEJUtQxhhjgpIlKGOMMUHJEpQxxpig1BgqSfxiItIbGIezQeJa4FpV/dlq52AmIgOB8aqa6n4dhbPx4wU4lTWeVtVHAxhirUTkZOAxoDNQCDyhqi+H2nUAiMiZwCNAe5xr+VuoXguAiCQBi4Gxqvp6KF6HiFwLvAzs82i+GXiX0LuWVsBLwIlAGfBPVf1LKP65VLEeVA3cf6j/A94HkoCHgYkikhjQwLwkIi4RuR6YCER5HHoAEKAjMAC4SkSuDECIRyUimcDHwF9x/gwuBR4VkVMJoeuAQ788PgLuUdUE4ELgGRHpR4hdi4dxQGuP16F4Hf2Ap1Q13uPxBqF5Lf/DKenWEhiME/NlhOa1AJaganMCEKmqz6jqAVV9D1gGXBzYsLz2AHATzi93T1cBD6vqDlXNBZ4Ebqzn2LzVDnhHVT9R1Qp37/U7YBihdR2o6haghap+JSJhQHPgILCLELsWABG5CkgElng0h9x1AP2BhUdoD6lrEZFBQAfgVlUtU9V1OL/DphBi1+LJhvhq1g1YUa0tB+gZgFh+iXGqOlZETqhqcA/JtAKWe7wvaK9JVX8Afqh6LSLNcAr//ocQuo4qqrpLROKAYpz/9x4HigixaxGR9sB9wFBggrstpP5uAYhIOM7w/RUi8jRQCryKM+QXUteCk2iXAPeLyNU4Q3wvAq8RetdyiCWomsXj/IX1VArEBSCWY6aqm4/QHO/+6XldIXFNItIU+AyYDcx3N4fcdeD84miC84vxS2Cvuz0krsX9S/0t4C5VzReRqkOh+HerBTAPeAM4D+iKM0xWNSQeStdS9Y+3qTg9qSycfzwUuY+H0rUcYgmqZnuA2GptccDuAMTiK3vcPz2vK+ivSUS64PziWA5czk/xh9R1AKhqBbAfmCci/wSy3YdC5Vr+Aqiq/rdae8j93VLVfOB4j6aFIvI8cLr7dchcC84kjxJVvd/9epGIvIozvAehdS2H2D2omi3HubHoKYvDu8ohRVV3APkcfl1BfU0iMgKn1/QpcIF7fD0Ur+N4EZlfrTkaCLVruQS4QER2ishOnKGiF3EmEYXSdSAi3UXkgWrNUTi93JC6Fpxhuzj35K4qEYTe36/DWA+qZlMAl4jcjjNF83ycYZlPAhpV3f0HuE9EFuMMy9wFPBvYkI5MRDoC44F7VfX5aodD5jrcFgKtReQOnDgHAdcB5+L8AgmJa1HVLM/XIrIQeMY9zXw3IXIdbjuBO0UkD+deTV/gVuAWnAlRoXQt3+AM5z0lInfiJKTrcCZKrSW0ruUQ60HVQFX343T1zwe2A/cC56hqUa0fDH5jgaU4/wPOxZnGPS6gEdXsZiABZ2r5bo/H44TWdaCqxcAZOPc6tgP/BK5X1amE2LXUIqSuQ1U3AWfhzGgrwYn3IVX9iNC7ljKc4coOOFPNJ+Css/uYELsWT7ajrjHGmKBkPShjjDFByRKUMcaYoGQJyhhjTFCyBGWMMSYoWYIyxhgTlCxBGWOMCUq2UNcYL4nI6/xUOuZIHsCptj4FSFDVeikn466PNx24UlVXHuH4LTi189r54Fyf42xP8V1dv8uYo7EelDHeuw2nMnQrnK0MAAZ6tD0JzHA/33OEz/vLrcCiIyUnP7gXeKlaSR1j/MJ6UMZ4yV0NohhARFLczUXuoqOeqr/2GxGJAf4IjKyP86nqYhHZAlwGvF4f5zSNlyUoY3zIvf/WoSE+EanE2Qn4jzj10eYBvwLuBq7AKbHzR1X9j/vzCcBTONtzVwKTgdtq2D4FnOKtO1V1qUcMA4DncWpHzsXp1XnGWLUjcW/3OWYBv1XVHBH5Etimqld4vH8sMEpVqyp//xenN/n6sf73MeZY2BCfMf73GPB7nG242wALcBLTAJxf9i+LSNV+Sv/ESWSn4tRWqwS+FpGa/jE5GvemgXBoU8cJwCKc7czfBG73ON4WZ1+tD3E25RyJs5fQE+63vAWcLSKe2zNcCrzt8XoC0Me9jb0xfmMJyhj/+4eqTlHVhTjV2XcDf1JVBZ7G2aunvYh0wOkRXaaqc929oiuAdsBpNXx3Nk4R0CoX4+wNdIuq5qjqaxze04nAmTDxlKquU9VZOJXhu7mPfwq4cBIfItIPpwDpRx7fsdZ9jv7H/F/CmGNgCcoY/1vt8bwUyFXVqirNZe6f0fyUJLSqcjuwDWcH3up7k1VpCWz1eN0DWKKqBzza5lQ9UdU1wH9F5B4ReVNE5gCPAOHu46U4vbpL3B+5DPhKVbd7fEcFzj5DqUe9cmPqwO5BGeN/B6q9rqjhfRHu9/bFGdrztP3nbz/0Xa5qbdVf7696IiI9gJk4+wdNBV7F2ZvqZo/3vwV85r4fdjFwxxHOGw6U1xCTMT5hPShjgscKIBJooqqrVXU1zt4+TwBdavhMPtDC4/VioLeIRHu09fN4fhXwo6qep6rPqur3QHsOT2qTcBLi3UAi8LnnCUUkDOe+Vb3NVjSNk/WgjAkSqqoi8hnwpojcjLND6sM4kytyavjYfJzZeFXeA+4DXhGRR3CS0/X81APbBmSJyHBgE3A2cIO7vSqOChF5FydBvefeDM9Td5yEtuCXXqsx3rAelDHB5Sqcqeif4kwRbwqcrKo7a3j/Fziz/YBDa7VOAtoCP+Js7/13j/c/B3yL0yuaD5wJ/AZIFZEMj/e9A8Rw+Oy9KiOA2Q1gd2kT5GxHXWNCmIjEAbnAaarqsx6NiIzGmfKe6Z4U4XlsJvBi1dotY/zFelDGhDD3rLsnOHySwy8mIu1F5ELgcWDcEZJTX5x7Xu/64nzG1MYSlDGh7+9ALxGpaSr6scgE/g3k4VS0qO4h4EZVPeiDcxlTKxviM8YYE5SsB2WMMSYoWYIyxhgTlCxBGWOMCUqWoIwxxgQlS1DGGGOC0v8Dfx3Ve7opTUEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "plot(results.r, label='r')\n",
    "\n",
    "decorate(xlabel='Time (day)',\n",
    "         ylabel='Distance from sun (million km)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
